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ABSTRACT 

The  Geophysical  Fluid  Dynamics  Laboratory  Modular  Ocean  Model  (GFDL 
MOM)  is  used  to  investigate  the  model  difference  between  compatible  and  incompatible 
surface  wind  and  buoyancy  forcing.  The  atmosphere  is  a  physical  system  in  which 
surface  wind  and  temperature  fields  are  related,  however  in  most  ocean  numerical 
models,  the  wind  stress  and  buoyancy  forcing  are  usually  specified  separately,  i.e.,  no 
constraint  between  the  surface  wind  stress  and  surface  air  temperature  is  considered.  In 
reality,  only  one  of  these  two  fields  can  be  prescribed  in  the  atmosphere-driven  ocean 
models.  When  the  surface  wind  field  is  prescribed,  the  surface  air  temperature  should 
be  derived,  and  vice  versa.  If  the  two  related  fields  are  treated  as  totally  independent  in 
forcing  the  ocean  models  the  results  will  be  distorted.  Since  the  model  solutions  depend 
upon  the  atmospheric  forcing,  it  is  important  that  we  study  the  compatibility  between  the 
wind  and  buoyancy  forcings  and  the  effect  which  incompatibility  might  have  on  the  ocean 
numerical  models. 

This  study  shows  that  the  surface  wind  and  buoyancy  forcing  widely  used  in 
ocean  numerical  models  are  incompatible.  Such  an  incompatibility  results  in  21  %  error 
in  the  total  northward  transport  of  heat,  16%  error  in  the  total  northward  transport  of 
salt,  25%  error  in  v  velocity,  and  16%  error  in  w  velocity. 
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I.     INTRODUCTION 

In  the  atmosphere  the  surface  wind  and  the  sea  surface  air  temperature  (SSAT)  are 
closely  related.  As  such,  thermal  field  and  the  surface  wind  field  cannot  be  prescribed 
independently  in  a  numerical  model. 

However,  in  current  ocean  numerical  models,  the  surface  wind  and  the  net  surface 
heat  flux  are  usually  prescribed  independently.  The  assumption  that  the  prescribed 
surface  wind  field  always  matches  the  surface  thermal  field  (or  the  net  surface  heat  flux) 
is  rather  unrealistic.  Contrary  to  the  assumption,  the  two  fields  do  not  match  each  other 
in  many  cases.  The  imbalance  between  the  surface  wind  and  the  surface  thermal  field 
gives  rise  to  certain  adjustment  processes,  thus  changing  the  original  forcing  field. 
Consequently,  during  the  entire  integration  of  the  ocean  numerical  models  these  two 
forcing  fields  will  not  maintain  the  prescribed  form. 

That  there  is  an  obvious  incompatibility  can  be  found  from  the  pioneering  work  by 
Takano  (1975),  which  states  that  the  SSAT  is  a  function  of  latitude  (y)  only,  and  that  the 
surface  wind  stress  is  a  function  of  position  (x,y)  and  season  (/) 

Ta  =  Ta(y)  ,       us  =  us(x,y,  t)  ,       vs  =  vs{x,y,  t)  (1) 

In  equation  (1)  the  atmospheric  surface  isotherms  are  given  by  straight  lines 
running  parallel  to  the  x-axis.  However,  surface  winds  have  both  x  and  y  components 
and  seasonal  variation.   It  is  thus  unlikely  that  the  atmospheric  surface  isotherms  would 


1 


remain  linear  and  in  a  steady-state  when  there  are  surface  winds  blowing  across  the 
isotherms. 

Haney  (1971)  states  that  by  employing  a  heat  budget  analysis  appropriate  to  zonally 
and  time  averaged  conditions  within  the  atmosphere,  it  can  be  shown  that  the  net 
downward  heat  flux  Q  at  the  ocean's  surface  is  expressed  as 


Q  =  cur.  -  Ta) 


(2) 


where  TA  is  an  apparent  atmospheric  equilibrium  temperature,  Ts  is  the  sea  surface 
temperature,  and  Q2  is  a  coefficient  determined  from  the  zonally  and  time  average  data. 
The  zonally  averaged  surface  zonal  winds  and  apparent  air  temperature  were  broadly 
used  as  mechanical  and  thermal  forcing  in  ocean  numerical  models  when  the  surface 
thermal  boundary  conditions  (2)  in  Haney  (1971)  are  met.  This  approach  is  well 
accepted  by  ocean  modelers. 

Several  combinations  of  surface  wind  stress  t  and  TA  in  ocean  models  are  listed  as 
follows  (Chu,  1992).  In  Davey  (1983)  model,  the  atmospheric  temperature  TA  is 
specified  by 


TAy)    = 


Ts   + 


(t„  -  rj 


1  -  cos 


« (y  -  ys) 


>    y  <  ys 

,     ys  <  y  <  ys  +  1 
,    ys  +  l  <  y        (3) 


with  /  =  2000  km,  Ts  =  12°C,  and  TN  =  8°C.   The  surface  winds  are  set  to  zero,  thus 


(us,    va)    =  0  (4) 

The  important  question  arises  as  to  whether  the  atmospheric  forcing  indicated  by 
(3)  and  (4)  can  last  2000  days,  the  period  over  which  the  model  was  integrated. 
Geostrophic  baroclinic  zonal  atmospheric  flow  is  usually  driven  by  the  north-south  SSAT 
gradient  while  the  cross  isothermal  surface  winds  is  generated  by  the  resulting  friction. 
Hence,  the  atmospheric  mechanical  forcing  given  by  (4)  is  not  consistent  with  the 
atmospheric  thermal  forcing  of  (3). 

There  are  other  renowned  ocean  model,  such  as  Robinson  et  al  (1977),  Semtner 
and  Mintz  (1977),  Cox  and  Bryan  (1984),  Cox  (1985),  Huang  (1989),  and  Willmott  and 
Darby  (1990).  In  general,  the  wind  stress  is  purely  zonal  (jy  =  0)  and  assumed  as  a 
sinusoidal  function  (or  other  form)  of  y.  TA  is  assumed  as  a  linear  function  of  y,  e.g., 
in  Robinson  el  al  (1977), 


JlliiL,       x     =  o,        TA  =  26    -  1 

yN-  ys       y  ^ 


xx  =  t0cos7c— ££,       t     =  0,       TA  =  26   -  -My  "  y8)  <5> 


Since  the  ocean  numerical  models  are  generally  integrated  over  long  time  interval, 
such  as  880  years  in  Bryan  and  Cox  (1984),  it  is  important  to  determine  whether  or  not 
the  surface  wind  and  the  SSAT  as  employed  by  the  ocean  numerical  models  are 
consistent.  Where  there  are  inconsistency,  the  effect  of  the  inconsistency  on  the  ocean 
circulation  should  be  further  investigated. 

The  main  purpose  of  this  study  is  to  address  these  differences.  For  this, 
experimentations  with  a  three-dimensional  primitive  equation  ocean  model  are  performed. 


Through  model  simulation  of  the  principal  aspects  of  the  Pacific  Ocean,  we  gain  a  better 
understanding  of  the  differences  between  the  atmospheric  parameters  representing 
mechanical  and  buoyancy  forcing,  i.e.,  the  surface  wind  and  SSAT  field. 

A  statement  of  the  problem  in  existing  ocean  numerical  models  is  contained  in  this 
chapter.  Chapter  II  gives  a  description  of  the  Modular  Ocean  Model  (MOM).  Chapter 
HI  describes  the  compatibility  between  surface  wind  stress  and  buoyancy  forcing. 
Chapter  IV  shows  the  errors  caused  by  incompatible  forcing.  Finally,  Chapter  V  states 
the  conclusions. 


H.    DESCRIPTION  OF  THE  MODULAR  OCEAN  MODEL  (MOM) 

A.      INTRODUCTION 

More  than  two  decades  have  passed  since  the  first  3-dimensional  primitive  equation 
numerical  ocean  model  was  used  to  study  the  most  basic  aspects  of  large-scale,  baroclinic 
ocean  circulation  (Bryan  and  Cox,  1967).  A  description  of  the  physics  and  numerics 
involved  was  published  in  Bryan  (1969).  In  this  model,  the  prediction  of  currents  is 
carried  out  using  the  Navier-Stokes  equations  with  three  basic  assumptions.  The 
Boussinesq  approximation  is  adopted,  in  which  density  differences  are  neglected  except 
in  the  buoyancy  term.  The  hydrostatic  assumption  is  made  in  which  local  acceleration 
and  other  terms  of  equal  order  are  eliminated  from  the  equation  of  vertical  motion. 
Lastly,  closure  is  attained  by  adopting  the  turbulent  viscosity  hypothesis  in  which  stresses 
exerted  by  scales  of  motion  too  small  to  be  resolved  by  the  grid  are  represented  as  an 
enhanced  molecular  mixing.  The  temperature  and  salinity  are  calculated  using 
conservation  equations,  again  utilizing  a  turbulent  mixing  hypothesis  for  closure.  The 
equations  are  linked  by  a  simplified  equation  of  state. 

For  the  purpose  of  computational  efficiency  several  techniques  are  used.  High 
speed,  external  gravity  waves  are  eliminated  by  the  "rigid-lid"  approximation,  and  a 
Laplacian  equation  is  solved  for  the  non-divergent,  vertically  averaged  flow.  The 
timestep  limitation,  i.e.,  the  half  pendulum  day  constraint  associated  with  inertia-gravity 
waves,  is  overcome  by  a  semi-implicit  treatment  of  the  Coriolis  term. 


Considerable  improvement  was  made  to  the  structure  of  the  FORTRAN  code  of 
this  model  by  Semtner  (1974)  who,  at  the  same  time,  added  various  features  to  the 
mathematical  formulation,  chief  of  which  was  the  use  of  "hole  relaxation"  (Takano, 
1974)  in  the  solution  of  the  external  mode  for  a  model  with  islands.  This  version  of  the 
model  has  been  adopted  by  many  investigators  and  has  seen  considerable  use  for  a 
number  of  years  in  the  ocean  modelling  work  at  GFDL.  During  this  time,  as  vector 
processing  machines  became  more  demanding  of  suitable  FORTRAN  structure, 
significant  changes  have  been  made  to  the  code  for  efficiency  purposes.  It  has  also  been 
generalized  in  several  ways,  among  which  is  the  incorporation  of  variable  grid  spacing 
in  the  horizontal,  and  an  arbitrary  number  of  tracer  prognostic  variables.  The  relaxation 
code  for  the  solution  of  the  external  mode  has  been  redesigned,  and  a  better  technique 
for  establishing  the  initial  guess  has  reduced  the  scans-to-convergence  considerably. 

Two  objectives  have  been  sought  in  designing  the  code  for  use  outside  of  GFDL. 
First,  it  has  been  made  universal  to  some  degree,  by  the  use  of  optional  code  lines.  A 
separate,  FORTRAN  coded  updating  utility  is  provided  to  carry  out  these  operations. 

B.       THE  MODEL  DESCRIPTION 

The  GFDL  Modular  Ocean  Model  (acronym  MOM)  version  1.0... Dec  1,  1990  is 
a  three  dimensional  primitive  equation  general  ocean  circulation  model  intended  for  use 
as  a  flexible  tool  for  exploring  ocean  and  coupled  air-sea  applications  over  a  wide  range 
of  space  and  time  scales. 


MOM  has  been  written  as  a  collaborative  effort  by  Ron  Pacanowski,  Keith  Dixon, 
and  Anthony  Rosati  at  the  National  Oceanographic  and  Atmospheric  Administration's 
Geophysical  Fluid  Dynamics  Laboratory  in  Princeton,  New  Jersey.  It  is  the  successor 
to  the  code  written  by  Michael  Cox,  documented  in  the  GFDL  Ocean  Tech  Report  #1, 
(1984).  As  was  the  case  for  the  Cox  model  and  the  Semtner  model  (UCLA  Dept.  of 
Meteorology  Tech.  Report  No.  9,  1974)  that  preceded  it,  MOM  is  a  Fortran 
implementation  of  equations  described  by  Kirk  Bryan  (1969). 

For  my  investigation,  I  use  the  Pacific  Ocean  as  the  region  of  interest  in  the  MOM 
model.  The  region  is  defined  by  160°  of  longitude  from  120°E  to  280°E  and  extending 
from  60°S  to  60°N.  A  constant  depth  of  5700  m  is  assumed.  The  horizontal  resolution 
is  5°  in  longitude  and  5°  in  latitude.  There  are  15  levels  in  the  vertical,  with  resolution 
varying  from  30  m  near  the  surface  to  836  m  near  the  bottom. 

The  time  manager  clock  parameters  for  setting  time  of  model  initial  conditions  are 
yearO  =  starting  year  (set  to  be  0),  monthO  =  starting  month  (set  to  be  1),  and  dayO  = 
starting  day  (set  to  be  1).  Using  logical  Julian  calendar  in  MOM  model,  the  time  period 
of  integrations  is  60  years.  The  time  step  for  both  barotropic  and  baroclinic  velocity  is 
1  hour.  For  density  and  tracers  (potential  temperature  and  salinity),  the  time  step  is  1 
day. 

The  MOM  model  uses  an  efficient  grid  system.  Horizontally,  tracer  quantities  are 
defined  at  the  centers  of  "t"  grid  boxes  and  velocities  are  defined  at  the  centers  of  "u,v" 
grid  boxes.  The  centers  of  "u,v"  grid  boxes  are  located  at  the  northeast  corners  of  "t" 
grid  boxes.    The  first  "t"  grid  box  is  located  in  the  southwest  comer  of  the  "t"  grid. 


This  grid  system  is  replicated  and  stacked  vertically  one  on  top  of  another  from  the 
surface  of  the  ocean  downward.  Vertically,  tracers  and  velocities  are  defined  at  the 
centers  of  their  corresponding  boxes  and  are  set  at  the  same  depths. 

The  values  for  the  mixing  coefficients  of  MOM  model  were  chosen  as  109  cm2/s 
and  1  cnr/s  for  the  horizontal  and  vertical  viscosity,  respectively,  and  2  x  107  cm2/s  and 
20  cnr/s  for  the  horizontal  and  vertical  diffusivity,  respectively.  The  bottom  drag 
coefficient  was  set  to  0. 

C.      BASIC  EQUATIONS  OF  MOM  MODEL 

The  basic  equations  of  the  model  are  written  here  in  continuous  form. 
Let 


m  =  sec<J> 
n  =  sin4> 
f  =  2Qsin<J> 


where  <j>  is  latitude. 


An  advective  operator, 


(6) 


T(\x)    =  m 


l, 


Wl+(v„i) 


w\l) 


(7) 


is  adopted  in  which  fi  is  any  scalar  quantity,  X  is  longitude  and  a  is  the  radius  of  the 

earth. 

The  equations  of  motion  are  then 


ut  +  r(u)    -  f v  =  -m—\-2-  |     +  Fu  (8) 

<3l  P0 


vt  +  r(v)    +  f  u  =  --[-^1     +  Fv  (9) 


where  p0  is  taken  to  be  unity. 

The  local  pressure,  p,  is  given  by  the  hydrostatic  relation, 


p(z)    =  ps  +   f  gpdz  (10) 

J  z 


where  ps  is  the  pressure  at  the  surface  of  the  ocean. 
The  continuity  equation  is 

ru)  =  o  (id 

The  conservation  equation, 

Tt   +  T(T)    =  FT  (12) 

applies  to  any  "tracer"  type  of  quantity  carried  in  the  model.    These  include  the  active 
tracers,  potential  temperature  and  salinity  (active  in  the  sense  that  they  appear  in  the 
equation  of  state),  and  any  passive  tracers  such  as  Carbon  14  or  Tritium. 
The  equation  of  state  is 

p   =  p(6/S,z)  (13) 


where  6  is  potential  temperature,  S  is  salinity  and  the  depth  dependence  arises  from 
compression  effects.   In  the  present  model,  (13)  is  represented  by  a  polynomial  fit  to  the 
Knudsen  formula  (Bryan  and  Cox,  1972). 
Let 


V2^    =  i7?2U-u    + 


m 


The  effects  of  turbulent  mixing  are 


(14) 


Fu  =  AMVuzz  +  A^a.-2  [V*u  +   (1  -  m2n2)  u  -  2nm2vk]  (15) 


Fv  =  Amvzz  +  A^-2  [W  +    (1   -  m2n2)  v  +  2nm2uk] 


(16) 


FT  = 


lTV 


L\ 


T. 


+  ATHa~2V2T 


(17) 


where  A  is  the  mixing  coefficient  corresponding  to  M  (momentum),  T  (tracer),   V 
(vertical),  and  H  (horizontal). 

In  nature,  vertical  mixing  is  known  to  be  a  rather  complex  function  of  vertical 
stability.  Since  this  process  is  still  not  well  understood,  we  have  adopted  a  simple, 
uniform  mixing  under  statically  stable  conditions,  and  an  infinite  mixing  under  statically 
unstable  conditions.    This  is  done  by  specifying  b  to  be 


6   =  \ 


*?    <0 
dz2 

*SL  >  o 

dz2 


(18) 
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At  lateral  walls,  the  boundary  conditions  are 


u,    v,    Tn  =  0 


(19) 


where  the  n  subscript  indicates  a  local  derivative  with  respect  to  the  coordinate  normal 
to  the  wall.    At  the  surface, 


pQAm(uz,    vz)    =    (t\    t*) 
Arv(Tz)    =  n 
w  =  0 


(20) 


The  "rigid-lid"  assumption  of  zero  vertical  motion  at  the  surface  filters  out  high  speed 
external  gravity  waves  which  would  otherwise  seriously  limit  the  length  of  the  time  step 
used  in  the  numerical  integration.  The  quantities  t\  t*  are  the  zonal  and  meridional 
components  of  the  specified  surface  stress,  and  r\  is  a  flux  through  the  surface,  of  the 
particular  tracer  involved. 
At  the  bottom, 


PoAMv(Uz>     Vzy>      =    ° 

r,  =  o 


w  =  -mua  1HX  -  va  1H^ 


(21) 


Combining  (8)  and  (9)  with  (10),  we  derive 


iut  =  ut'  -  ma  ~xpx 


Wt  =  vt'  -  a~xp£ 


(22) 
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where 


uj  =  -F(u)    +  f v  -  mga'lf°pkdz'  +  Fu  (23) 


vj  =  -T(v)    -  fu  -  a'1f°p(^dz/  +  Fv  (24) 


Let  us  define 


u  =  u  +  u;  v  =  v  +  v  (25) 


where 


]I  =  H^C^dz  (26) 


Then,  since  /?*  is  not  a  function  of  depth, 

at  =  ut   ~  «?;  ^t  =  vt  -v?  (27) 

Since  all  terms  on  the  right  of  (23)  and  (24)  are  known,  (27)  may  be  solved  for  the 
internal  modes  of  momentum.  Under  the  rigid-lid  boundary  condition,  the  external  mode 
of  momentum  may  be  represented  by  a  volume  transport  stream  function,  \p, 

u  =  -(afl)"1^;  v  =  m(aH)  "1i|rA  (28) 

This  is  shown  by  integrating  (11)  vertically,  substituting  (28)  and  noting  that  the 
boundary  conditions  (20)  and  (21)  on  w  are  satisfied.  A  prognostic  equation  for  \p  may 
be  obtained  by  averaging  (22)  vertically,  and  eliminating  terms  in/f  by  applying  the  curl 
operator, 
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Substituting  (28), 


curlz(vt,    ut)    =  ma 


-i 


'»)». 


(29) 


/m|r 


rX 


Ha 


* 


t<t> 


mHa 


=  a 


J<t) 


v' 


tA 


V   m 


(30) 


The  boundary  condition  on  4>  at  lateral  walls,  corresponding  to  (19)  is 


♦♦  =  *h  =  ° 


(31) 


This  condition  is  satisfied  by  setting  $  constant  over  each  unconnected  land  mass 
comprising  the  ocean  boundary.  In  the  case  of  an  enclosed  basin  with  no  islands,  yp  may 
be  set  to  zero  over  the  boundary  forming  land  mass.  If,  in  addition,  islands  are  present, 
the  associated  constant  for  each  island  reflects  the  net  flow  around  the  island  and  must 
therefore  be  predicted  by  the  governing  equations.  The  method  used  is  "hole  relaxation" 
in  which  the  line  integral  of  the  quantity  Vps,  taken  around  the  island,  is  required  to 
vanish.  Averaging  (22)  vertically,  integrating  around  the  coast  of  the  island  and  setting 
the  contribution  due  to  ps  to  zero,  the  predictive  equation, 


f 


mty 


a 


H 


d<J)  - 


* 


t* 


mH 


dX 


£  v'  cd$ 


u' 


m 


d\ 


(32) 


is  obtained.    Applying  the  Stokes  theorem  yields  a  more  useful  form, 


1 


MVtk 

H 

+ 

■dA  =  a 


^  -  (5), 


dA 


(33) 


Note  that  (33)  is  simply  an  area  integral  of  (30),  taken  over  the  island. 
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D.      PROGRAM  FLOW  OF  MOM  MODEL 

The  FORTRAN  code  of  the  model  consists  of  the  main  program  OCEAN,  and 
seven  subroutines.    Their  functions  are  described  below. 

OCEAN:  Performs  all  operations.  This  is  done  only  once  at  the  beginning  of 
each  run  of  the  model.  The  program  calls  STEP  once  per  timestep, 
and  attends  to  operations  which  must  be  done  at  the  end  of  each  run. 

STEP:  Called  once  per  timestep  by  OCEAN,  it  initializes  various  quantities, 

bootstraps  the  row-by-row  computation  of  prognostic  variables, 
manages  the  I/O  for  the  latter,  and  performs  various  analysis 
procedures  on  the  progressing  solution. 

CLINIC:  Called  once  per  row  by  STEP,  it  computes  the  internal  mode 
component  of  the  u  and  v  velocities  as  well  as  the  vorticity  driving 
function  for  use  by  RELAX  later  in  determining  the  external  mode. 

TRACER:  Called  once  per  row  by  STEP,  it  computes  temperature,  salinity,  and 
any  tracers  which  are  carried  in  the  model. 

RELAX:  Called  once  at  the  end  of  each  timestep  by  STEP,  it  takes  the  vorticity 
driving  function  computed  in  CLINIC  and,  using  sequential 
overralaxation,  solves  the  Laplacian  equation  for  the  external  mode  of 
velocity  in  terms  of  a  mass  transport  stream  function. 

STATE:  Called  by  CLINIC  and  TRACER  once  per  row,  and  STEP  in  the 
bootstrap  procedure,  it  computes  normalized  densities  by  using  a  third 
order  polynomial  fit  to  the  Knudsen  formula. 
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MATRIX:     Called  by  STEP  on  specified  timesteps,  it  is  a  general  2-dimensional 

array  printing  routine. 
ODAM:        (Ocean  Direct  Access  Manager) 

Used  primarily  by  STEP,  it  is  a  set  of  routines  which  is  responsible 

for  handling  the  transfer  of  data  between  memory  and  disc  (virtual 

disc  residing  in  memory  in  the  core  contained  mode). 

In  summary,  for  one  timestep,  OCEAN  calls  STEP,  which  establishes  the  proper 

data  in  memory  and  calls  CLINIC  and  TRACER  row-by-row  from  south  to  north  through 

the  basin.    Upon  completing  the  fmal  row,  STEP  calls  RELAX  and  returns  control  to 

OCEAN,  which  may  call  STEP  for  another  timestep. 

Two  additional  subroutines,  FINDEX  and  FILTER  are  used  when  Fourier  filtering 
is  needed  to  overcome  the  timestep  limitation  arising  from  convergence  of  meridians  at 
high  latitudes. 

E.      DISC  I/O  SYSTEM 

Next  to  the  equations  themselves,  the  component  of  the  code  which  accounts  for 
the  greatest  complexity  is  the  I/O  system.  The  purpose  of  this  system  is  twofold.  First, 
having  a  complete  record  of  all  prognostic  variables  on  permanent  disc  at  the  end  of  each 
timestep  allows  restarting  an  experiment  from  an  earlier  run,  or  from  a  machine 
malfunction.  Secondly,  during  a  run,  it  is  generally  impossible  to  fit  the  three 
dimensional  arrays  in  memory  entirely.  The  I/O  system  is  designed  to  feed  data  to  and 
from  memory  in  a  row-by-row  manner.    Data  for  one  row,  including  all  east- west  and 
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vertical  grid  points  is  termed  a  "slab".  At  any  one  time,  only  the  slabs  necessary  for  the 
computation  are  present.  While  the  computation  for  one  row  proceeds,  the  I/O  system 
is  feeding  the  slab  just  computed,  back  to  disc  and  fetching  the  slab  which  will  be  needed 
to  compute  the  next  row.  If  the  disc  transfers  are  done  fast  enough,  no  wait  for  data  will 
be  needed  at  the  beginning  of  each  row  and  the  system  is  said  to  be  completely 
"buffered". 

Three  disc  units  are  needed  for  the  process  when  using  centered  differencing  in 
time;  one  for  the  N  timestep  data  (read),  one  for  the  N-l  timestep  data  (read),  and  one 
for  the  newly  computed  N+ 1  timestep  data  (written).  Since  the  N  level  on  one  timestep 
becomes  the  N-l  level  for  the  next,  it  is  convenient  to  permute  the  disc  units  to  minimize 
data  transfers.  Thus,  on  timestep  1,  N-l  and  N  data  are  read  from  units  13  and  14,  and 
the  N+l  data  is  written  to  15.  On  timestep  2,  units  14  and  15  are  read  and  the  N+l 
data  is  written  to  13.  On  timestep  3,15  and  13  are  read  with  the  N+ 1  data  going  to  14, 
etc. 

Possibly  the  most  abstruse  feature  of  the  model  is  the  manner  in  which  the  "slab 
incidental  data"  is  handled.  Just  as  the  slab  system  reduces  the  row  dimension  of  the  3- 
dimensional  prognostic  variables  to  3,  it  can  also  be  used  to  reduce  the  row  dimension 
of  2-dimensional  variables  which  would  otherwise  add  considerably  to  the  memory 
requirements.  Also,  if  these  variables  are  constant  in  time,  there  is  no  need  to  keep 
multiple  time  records  of  them.  Consider  two  arrays,  A  and  B,  for  which  there  is  data, 
invariant  in  time,  at  each  row  and  column  horizontally  across  the  grid.  They  may  be 
carried  in  the  model  as  "slab  incidental  data",  thereby  reducing  their  row  dimension  to 
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3.  Furthermore,  let  A  reside  in  the  slab  corresponding  to  the  N-l  time  level  of  the 
primary  slab  data,  and  B  reside  in  the  N  level  slab.  The  diagram  below  illustrates  the 
6  timestep  cycle  which  the  permuting  disc  units  execute,  with  N-l,  N  and  N+ 1  denoting 
the  primary  slab  data,  and  A,  B  denoting  slab  incidental  data. 

TIMESTEP:  12  3  4  5  6  7 

UNIT 

13  N-l,  A     N+1,B      N,B      N-1,B     N  +  1,A      N,A      N-1,A 

14  N,B      N-1,B     N+1,A      N,A      N-l, A     N+1,B      N,B 

15  N+l.A      N,A      N-l, A     N+1,B      N,B      N-1,B     N+1,A 

Note  that  on  even  timesteps,  the  slab  incidental  data  enters  memory  in  the  wrong 
slabs  and  must  be  switched  between  slabs  N  and  N-l.  Also,  storing  of  A  and  B  into  the 
N+l,  slab  must  be  alternated  on  successive  timesteps.  In  the  base  code,  FKMU  and 
WSY  are  "A"  type  arrays,  and  FKMT  and  WSX  are  "B"  type  arrays. 

Disc  unit  12  is  used  primarily  for  the  storage  of  2 -dimensional,  horizontal  fields. 
It  is  divided  into  7  +  blocks  in  the  following  manner: 

Blocks  1-3:  permuting  blocks  for  the  stream  function  at  N-l ,  N,  and  N+ 1  time 
levels 

Block  4:  reciprocal  of  total  depth 

Blocks  5-6:       permuting  blocks  for  former  relaxation  solutions 
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Block  7  +  :  start  and  end  indices  (additional  blocks  are  added  as  necessary 
when  filtering  indices  fill  block  7) 

Finally,  unit  11  contains  the  timestep  counter,  total  elapsed  time,  and  the  area  and 
volume  of  the  basin. 

Data  is  fed  to  and  from  these  units  by  means  of  the  entry  points  of  subroutine 
ODAM.  It,  in  turn,  must  use  a  direct  access  I/O  package  provided  locally.  If 
FORTRAN  direct  access  I/O  is  available,  with  a  facility  to  buffer  the  operations  (such 
as  a  FIND  statement),  the  QDAM  calls  in  ODAM  of  the  base  code  may  be  replaced  by 
the  appropriate  FORTRAN  statements.  Otherwise,  a  specially  written  set  of  I/O  utilities 
such  as  QDAM  must  be  supplied.  If,  instead,  the  core-contained  option  is  invoked,  no 
such  package  in  needed. 
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HI.  COMPATIBILITY  BETWEEN  SURFACE  WIND  STRESS  AND 
BUOYANCY  FORCING 


Wind  stress  and  buoyancy  forcing  are  independently  designated  in  most  ocean 
models,  e.g.,  no  constraint  has  been  considered  between  the  surface  wind  stress  and 
surface  air  temperature  (or  net  surface  heat  flux).  The  atmosphere  actually  is  a  unified 
physical  system.  For  atmospheric  motion  on  the  scale  of  the  ocean  basin  considered 
here,  the  wind  and  temperature  fields  are  related  to  each  other  through  the  thermal  wind 
relation.  Only  one  of  the  two  (either  wind  or  air  temperature)  can  be  prescribed  in  the 
atmosphere-driven  ocean  models.  Thus,  treating  the  2  fields  as  totally  independent  will 
produce  physically  unrealistic  atmospheric  forcing  in  the  ocean  models.  This  gives  rise 
to  a  distortion  of  model  results  since  the  solutions  largely  depend  upon  the  atmospheric 
forcing.  Therefore,  it  is  important  to  study  the  compatibility  between  the  wind  and 
buoyancy  forcing  and  the  effect  of  incompatibility  on  the  model  results  before  running 
any  ocean  numerical  model. 

In  many  ocean  models  (e.g.,  Davis,  1983,  Robinson  et  al.  1977,  etc.)  the  SSAT 
is  assumed  zonal  symmetric.  Therefore,  if  the  surface  wind  («,*,  vs*)  and  SSAT  fields 
used  in  ocean  models  are  compatible  then  they  can  be  denoted  by  combinations  of 
sinusoidal  functions  (Chu,  1992),  i.e., 


Us*    =    UtY,  Usn   Sin   kn*y  (34) 


19 


vs*   =   c/rE  vsn  sin  K*V  (35) 

n 

6S.  =    (AD    "£Qsn  cos  ^n7ty  (36) 

where  the  subscript  '*'  indicates  the  dimensional  variables,  kn  (n  =  1,  2,...)  are  not 
necessary  integers.    It  should  be  noted  that  («/,  vs*)  is  the  total  surface  wind,  i.e., 

us.  =  ug*  +  u*>       vs*  =  V*>       at  z  =  0  (37) 

Since  the  flow  variables  are  independent  of  x  (zonally  symmetric),  the  geostrophic 
flow  should  be  purely  zonal,  therefore,  at  the  APBL  top  the  geostrophic  wind  would  be 

uG<   =  uGm(y)  ,       vG,  =  0  (38) 

The  scale  of  thermally  forced  surface  winds  (Chu,  1989)  is 

UTS-g£T  (39) 

T        2QLy60 

where  b  =  (v/0)'A,  is  the  Ekman  depth,  v  is  the  atmospheric  vertical  eddy  viscosity,  fl 
is  the  angular  velocity  of  the  earth's  rotation,  and  AT/Ly  is  the  characteristic  SSAT 
gradient.  If  5  =  1  km  (Holton,  1972),  AT  =  14°k,  Ly  =  2000  km  (Robinson  et  al, 
1977),  and  80  =  288 °k,  the  scale  for  the  thermally  driven  surface  wind  (UT)  approaches 
3  m/s,  and  the  scale  of  the  geostrophic  wind  at  the  APBL  top  (£/)  can  reasonably  be 
assumed  to  be  10  m/s.  Therefore,  if  we  completely  neglect  this  thermally  driven  surface 
wind  in  ocean  models,  we  will  make  almost  30%  error  in  the  surface  wind  forcing. 
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Based  on  the  K-closure  APBL  model  proposed  by  Kuo  (1973),  Chu  (1992)  found 
these  wind  and  temperature  fields  are: 

u{y,z)    =  Y^  [-bn  exp(-v/2Y  z) 

n 

+  2fo  t      f"^'    exp(X„jZ)  ]    sin  knny 

j  -  1     Anj  ZA 


6(y,z)    =  T  [^^"expC-v^y  z) 

3 


+  2knnRi  ]£       2    "J expUnj.z)]    cos  kntzy  (41) 

j  =  i   Xnj  -  2v 


v(y,z)    =  -||   =  -]T  £   a^A^exptt^z)    sin  jcany  (42) 


n    j  =  1 


evaluated  at  z=0,  i.e., 


"sn    =     -*n    +    2^0   £     TT^     +    VUGn.  *    =    ^     2,    .    .    .  (43) 


v„--S  a»A;  (44) 

J  =  1 


es„  ■  ^^  +  2k**Si  t  T^~  <45> 


*-»  M    *£,   ■■   2Y 


\n 


where  7  is  the  nondimensional  decay  rate  (7  =  0.01,  denoting  100  day  decay  time  scale) 
for  Rayleigh  friction  and  Newtonian  cooling,  f0  =  sin<£,  is  the  nondimensional  Coriolis 
parameter,  taken  as  locally  constant  except  very  close  to  the  equator,  <f>  is  latitude,  and 
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Ri   =  Ri{b/Ha)2  (46) 

is  the  generalized  Richardson  number.    Ha  =  10  km,  is  the  thickness,  and 

H    N 

Risiii^)2  (47) 

is  the  Richardson  number.   Here,  N  is  the  Brunt- Vaisala  frequency.   A  ■(/'  =  1,2,3)  are 
the  eigenvalues.    anl,  an2,  and  a^  are  integral  constants.    \p  is  the  streamfunction. 

The  boundary  condition  in  the  vertical  are  listed  as  follow.     The  dependent 
variables  should  remain  finite  as  z  -»  oo,  i.e., 


lim(  |i|/|, 


dty 


,  \u\,  lei)  <  ~  (48) 


dz 
The  surface  boundary  conditions  (at  z  =  0)  are 


*|z=0=0,      VLUa+u=Mf±l      v  =  M^z,      0--M1|§=0s         (49) 


where  M,  M1  are  given  constants,  and 

\x  =  U/UT  =3.33  (50) 

the  parameter  M  =  0.5,  which  is  a  measure  of  the  effective  depth  of  the  constant  stress 
sub-layer  (Kuo,  1973). 

The  surface  boundary  condition  (48)  and  (49)  leads  to  (assuming  Mr  =  0) 


t   *nj  •  °  <51) 

J     =    1 
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Usn   =  j2yMbn   +  2f0MY; 


Anj    ~    2Y 


(52) 


£     ^(1    "    ^nj^nj    =    0,       01        £    A 


nj  ™-nj 


J  =  1 


J  =  1 


Vsn 
M 


(53) 


Solving  the  linear  algebraic  equation  (44),  (51),  and  (53)  with  respect  to  anl,  an2,  and  an3, 
we  obtain 


where 


I 


nl 


lnl 


V. 
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sn '  n2 
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1         sn 
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(56) 


(57) 


h    = 


1 

1 

1 

Ki 

An2 

K3 

X2 

Anl 

X2 

An2 

X2 

An3 

Elimination  of  bn  from  (43),  (45),  and  (52)  leads  to 


(58) 
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where 


where 


Usn    =    -RlnQsn    +    #l^Gn 


(59) 


Vsn    =    ~R2nQsn    +   H2nUGn 


(60) 


kjK, 


^2n 


5 (1   +  ^2y  M)  — 


2Y   ^ 


JV 


(61) 


k„izM 


Rln    ~    InR2n 


(62) 


(63) 


^  s 


2Mkln2Ri 


4MfQy)   £ 


^ 


f?x  Jn(Xi,-  -  2Y: 


(64) 


*n    S 


2kl%2Ri  , 

( 1    +  v/2y  M)    -  4 f 0itfy 


^ov^Y 


2^oE 


"nj 


/ra    ln(X^  -  2y) 


1nj^nj 


lna2nj-2y)  (65) 

The  solutions  (59),  and  (60)  show  that  the  surface  winds  (usn,  vsn)  are  driven 

mechanically  by  the  geostrophic  winds,  uGn,  and  thermally  by  the  surface  thermal 

conditions,  -Rln0sn,  and  -R2Jsn. 

The  coefficient  in  (59)  and  (60)  Rln,  R2n,  HIn,  and  H^  for  different  kn  can  be 

computed  by  using  the  dependence  of  eigenvalues  Xn/,  Xn2,  and  X„j  on  kn  for 
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Ri   =  0.01  (66) 

Non-zero  values  for  these  coefficients  indicate  that  only  two  can  be  prescribed 
among  the  four  variables  usn,  vsn,  Bsn,  and  uGn,  the  other  two  should  be  derived  from  (59) 
and  (60).  Thus,  (59)  and  (60)  can  be  treated  as  a  compatibility  condition  among  these 
four  variables.  Obviously,  it  is  over- specified  if  we  prescribe  both  surface  wind  stress 
and  buoyancy  forcing  in  ocean  numerical  models. 

We  can  use  (59)  and  (60)  to  verify  the  consistency  between  buoyancy  forcing  (3) 
and  surface  wind  (4)  which  is  used  in  the  Davey  (1983)  model.  Here,  TA(y)  j±  0,  the 
Fourier  cosine  series  of  TA  is 

rA(y)   =  (TN  -  r5)£esn  cos  imy  (67) 

n 

The  surface  winds  (us,  vs)  are  set  to  zero,  therefore,  the  components  of  their  Fourier  sine 
series  should  also  be  zero,  i.e., 


"sn    =    0'  vsn    =    0  (68) 


Non-zero  values  of 


B     (HlnR2n  -  H2nRln) 
H2n 

shows  that  the  linear  algebraic  equations  (59)  and  (60)  about  6sn  and  uGn  only  have  zero 
solutions,  i.e.,  $sn  =  0,  uCn  =  0,  when  usn  =  0,  vsn  =  0.  Therefore,  the  thermal  forcing 
(3)  is  not  compatible  with  the  surface  wind  forcing  (4). 
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Purely  zonal  surface  winds  are  commonly  used  in  ocean  models,  such  as  the 
surface  boundary  conditions  (5).  Here,  three  variables  us,  vs,  and  TA  are  prescribed.  We 
expand  these  surface  variables  into  Fourier  sine  and  cosine  series 


(",*    vs)    = 


x  -r^E  <"«'    vsn)    sin  rniy,  vsn  =  0           <™> 

ta  =    (AD^esn  cos  rrny  (71) 

n 

ug  =  UY,  uGn  sin  nny  (72) 


n 


The  Fourier  components  usn,  vsn,  6m,  and  uGn  should  satisfy  the  two  linear  algebraic 
equations  (59)  and  (60): 


USn    =    "*lAn    +    HlnUGn  (73) 


0    =    -RzrPsn   +  Hzn*Gn  <74> 


which  shows  that  purely  zonal  surface  winds  appear  only  when  the  surface  thermal 
forcing  is  balanced  by  the  geostrophic  forcing 

*sn    =    ^UGn  (75) 

Rsn 

Elimination  of  uGn  from  these  two  equations  leads  to  a  relationship  between  usn  and 
$sn,  (Chu,  1992): 

u      =-(0  0  (76) 

sn  th    sn 
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which  is  the  compatibility  condition  for  zonally  symmetric  zonal  wind  and  SSAT  fields. 
Independently  prescribing  surface  wind  and  SSAT  fields,  such  as  in  (5),  violates  this 
condition  and  therefore,  is  inconsistent. 
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IV.    ERRORS  IN  NUMERICAL  INTEGRATION  OF  THE  MOM  MODEL 
CAUSED  BY  INCOMPATIBLE  FORCING 


A.      SURFACE  BOUNDARY  CONDITIONS 

The  MOM  model  and  its  parameters  have  been  discussed  in  Chapter  n.  Two  kinds 
(compatible  and  incompatible)  of  surface  boundary  conditions  (SST,  salinity,  wind  stress 
in  x  direction,  and  wind  stress  in  y  direction)  were  used  in  MOM  model  to  find  the 
differences  between  the  results  of  the  two  conditions. 

The  incompatible  surface  boundary  condition  (SST,  salinity,  wind  stress  in  x 
direction,  and  wind  stress  in  y  direction)  which  is  used  in  the  general  circulation  model 
by  Marotzke  and  Willebrand  (1991)  is  shown  in  figure  la.  The  zonal  (eastward)  wind 
stress  (solid  curve)  and  the  apparent  atmospheric  temperature  in  degrees  Celsius  (dashed 
curve)  are  functions  of  latitude  and  symmetric  about  the  equator,  so  they  are  only  shown 
for  one  hemisphere. 

The  sea  surface  temperature  (SST)  in  °C  is  symmetric  about  the  equator,  and 
follows  a  cosine  law  in  latitude  with,  a  difference  of  27 °C  between  high  and  low  latitudes 
(Fig.  lb).  The  salinity  is  analogous  to  the  sea  surface  temperature  and  symmetric  about 
the  equator.  The  salinity  also  follows  a  cosine  law  in  latitude  with  amplitude  2.5  psu 
(Fig.lc). 
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Many  independently  and  simple  prescribed  wind  and  buoyancy  forcing  in  ocean 
models  have  a  similar  form  as  used  by  Marotzke  and  Willebrand  (1991)  in  a  steady-state 
wind  and  thermally  forced  ocean  circulation  model  (Fig.  1,  2  and  3). 

The  latitudinal  surface  atmospheric  temperature  gradient  generates  extra  surface 
winds  (Mp  vr).  The  scale  of  this  thermally  driven  surface  winds  is  given  by  (39). 
Corresponding  to  SSAT  field  described  in  figure  lb,  the  thermally  driven  surface  winds 
can  be  predicted  by  a  K-closure  Atmospheric  Planetary  Boundary  Layer  (APBL)  model 
(Chu,  1989).    From  wind  and  temperature  fields  (40),  (41)  and  (42),  we  obtain 


uT 


vT 


UT 

Grp 


cosa    -since 
since    cosa 


ve  I  <77> 


where  a  is  the  deflection  angle  (angle  between  VT  and  V0a  |  z=0).  The  wind  field  with 
the  thermally  driven  component  (i.e.,  u+uT,  vT)  is  illustrated  by  solid  curves  and  without 
the  thermally  driven  component  is  illustrated  by  dashed  curves  in  figure  2  and  3. 

From  figure  2  and  3,  we  compute  the  Ekman  current  change  due  to  compatible 
forcing  (Fig.  4).  The  result  in  a  stronger  westward  and  poleward  current  between  about 
5°  and  30°  latitude  in  both  hemisphere.    We  also  compute  the  wind  stress  curl  for  the 
wind  field  with  and  without  thermally  driven  component. 
Wind  stress  curl  for  forcing  without  thermally  driven  component: 

K  •  curl  f  =  -^   -  ^  (78) 

ox         ay 

t  is  independent  of  x,  so 
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dz 
K  •  Curl  x  =  --— 
dy 


(79) 


Wind  stress  curl  for  forcing  with  thermally  driven  component: 


K  ■  Curl  x   =  -^   -  ^ 
ox         dy 


(80) 


J*,  is  independent  with  x,  so 


K  •  Curl  T  =  - 


3** 

dy 


(81) 


The  difference  between  the  wind  stress  curl  for  wind  field  with  and  without  thermally 
driven  component  is  not  negligible. 

B.      ERRORS    CAUSED    BY    INCOMPATIBLE    WIND    AND    BUOYANCY 
FORCING 

We  consider  a  rectangular  ocean  of  uniform  depth  in  the  MOM  model.  The 
thickness  (m)  and  depth  (m)  in  the  vertical  of  the  15  level  MOM  model  is  specified  as 
below: 


Level 

Thickness  (m) 

Depth  (m) 

1 

30.00 

30.00 

2 

46.15 

76.15 

3 

68.93 

145.08 

4 

99.93 

254.01 

30 


5 

140.63 

385.64 

6 

192.11 

577.75 

7 

254.76 

832.51 

8 

327.95 

1160.46 

9 

409.81 

1570.27 

10 

497. 1 1 

2067.38 

11 

585.36 

2652.74 

12 

669.09 

3321.83 

13 

742.41 

4064.24 

14 

799.65 

4863.89 

15 

836.10 

5699.99 

1.      The  total  northward  transport  of  heat 

The  total  northward  transport  of  heat,  which  was  calculated  by  compatible 
surface  wind  stress  and  buoyancy  forcing,  after  3,  6,  9,  and  12  years  is  shown  in  figure 
5.  The  difference  in  the  total  northward  transport  of  heat  between  compatible  and 
incompatible  surface  wind  stress  and  buoyancy  forcing  after  3,  6,  9,  and  12  years  is 
provided  in  figure  6.  In  all  plots  the  difference  is  computed  from  the  compatible  forcing 
run  minus  the  incompatible  forcing  run. 

The  total  northward  transport  of  heat,  which  was  calculated  by  compatible 
surface  wind  stress  and  buoyancy  forcing,  after  51,  54,  57,  and  60  years  is  shown  in 
figure  7.  The  difference  in  the  total  northward  transport  of  heat  between  compatible  and 
incompatible  surface  wind  stress  and  buoyancy  forcing  after  51,  54,  57,  and  60  years  is 
provided  in  figure  8. 
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When  the  northward  transport  of  heat  by  compatibility  consideration  is  large 
in  figure  5,  the  difference  in  the  northward  transport  of  heat,  shown  in  figure  6,  due  to 
the  incompatibility  forcing,  is  also  large.  The  error  reaches  to  21%  after  51  year's 
integration.  The  difference  in  the  northward  transport  of  heat  between  the  compatibility 
and  incompatibility  consideration  increases  with  the  integration  time. 

In  northern  (southern)  hemisphere,  an  additional  northeasterly  (southeastly) 
wind  stress  exists  in  the  compatible  forcing  case  (Fig.  2,  3).  Therefore,  an  extra 
northwest  (southwest)  Ekman  transport  should  appear  when  we  use  the  compatible 
forcing  (Fig.  4). 

2.      The  total  northward  transport  of  salt 

The  total  northward  transport  of  salt  which  was  calculated  by  compatible 
surface  wind  stress  and  buoyancy  forcing,  after  3,  6,  9,  and  12  years  is  shown  in  figure 
9.  The  difference  in  the  total  northward  transport  of  salt  between  compatible  and 
incompatible  surface  wind  stress  and  buoyancy  forcing  after  3,  6,  9,  and  12  years  is 
provided  in  figure  10. 

The  total  northward  transport  of  salt  which  was  calculated  by  compatible 
surface  wind  stress  and  buoyancy  forcing,  after  51,  54,  57,  and  60  years  is  shown  in 
figure  1 1 .  The  difference  in  the  total  northward  transport  of  salt  between  compatible  and 
incompatible  surface  wind  stress  and  buoyancy  forcing  after  51,  54,  57,  and  60  years  is 
provided  in  figure  12. 

When  northward  transport  of  salt  by  compatibility  consideration  is  large  in 
figure  9,  the  difference  in  the  northward  transport  of  salt  in  figure  10,  due  to  the 
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incompatibility  forcing,  is  also  large.  The  error  reaches  to  16%  after  51  year's 
integration.  The  difference  in  the  northward  transport  of  salt  between  the  compatibility 
and  incompatibility  consideration  increases  with  the  integration  time. 

3.      The  meridional  mass  transport 

The  meridional  mass  transport,  which  was  calculated  by  compatible  surface 
wind  stress  and  buoyancy  forcing,  after  30  years  is  shown  in  figure  13.  The  difference 
in  the  meridional  mass  transport  between  compatible  and  incompatible  surface  wind  stress 
and  buoyancy  forcing  after  30  years  is  provided  in  figure  14. 

The  meridional  mass  transport,  which  was  calculated  by  compatible  surface 
wind  stress  and  buoyancy  forcing,  after  60  years  is  shown  in  figure  15.  The  difference 
in  the  meridional  mass  transport  between  compatible  and  incompatible  surface  wind  stress 
and  buoyancy  forcing  after  60  years  is  provided  in  figure  16. 

From  above,  the  latitudinal  distribution  of  meridional  mass  transport  after  60 
years  integration  indicates  that  the  compatibility  consideration  meridional  mass  transport 
is  larger  than  that  of  incompatibility  consideration  in  all  latitudes  of  the  northern 
hemisphere  and  smaller  in  all  latitudes  of  southern  hemisphere.  The  difference  between 
the  two  consideration  (compatibility  and  incompatibility)  increases  with  the  latitude. 

In  figure  16,  the  difference  of  mass  transport  is  strong  between  25 °S  and 
25  °N  due  to  the  extra  Ekman  transport  when  the  compatible  forcing  is  used  (Fig.  4). 
The  error  is  more  then  5  Sverdrup  in  northward  mass  transport.  In  deeper  sea  water  it 
is  about  5  Sverdrup  less. 
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4.  Meridional  section  profile  of  0 

The  meridional  section  profile  of  0  at  160°E  after  60  years  is  shown  in  figure 
17  and  the  difference  is  provided  in  figure  18.  The  meridional  section  profile  of  0  at 
160°W  after  60  years  is  shown  in  figure  19  and  the  difference  is  provided  in  figure  20. 
The  meridional  section  profile  of  $  at  120°W  after  60  years  is  shown  in  figure  21  and 
the  difference  is  provided  in  figure  22. 

From  the  difference  of  meridional  section  profile  of  0  in  figure  18,  20,  and 
22,  we  can  see  that  the  potential  temperature  0  is  always  lower  in  the  compatible  forcing 
than  in  incompatible  forcing  near  the  equator.  The  errors  caused  by  incompatible  forcing 
in  upper  is  larger  than  that  in  the  deeper  sea  water. 

The  mass  transport  is  northward  in  north  hemisphere  and  southward  in  south 
hemisphere.  Thus  divergence  occurs  near  the  equator.  So  the  cold  deeper  sea  water 
moves  upward  toward  the  surface  to  make  the  temperature  lower  near  the  equator. 

5.  Meridional  section  profile  of  v  velocity 

The  meridional  section  profile  of  v  velocity  at  160°E  after  60  years  is  shown 
in  figure  25,  and  the  difference  is  provided  in  figure  26.  The  meridional  section  profile 
of  v  velocity  at  160°W  after  60  years  is  shown  in  figure  27,  and  the  difference  is 
provided  in  figure  28.  The  meridional  section  profile  of  v  velocity  at  120°W  after  60 
years  is  shown  in  figure  29,  and  the  difference  is  provided  in  figure  30. 

From  the  meridional  section  profiles  of  v  velocity  in  figure  25,  27,  and  29, 
the  v  velocity  is  positive  from  the  equator  to  about  25  °N  where  the  velocity  changes 
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direction  (opposite  sign)  toward  high  latitude.  The  v  velocity  in  the  southern-hemisphere 
is  compared  with  the  relation  in  the  northern-hemisphere. 

From  figure  14  and  16,  the  latitudinal  distribution  of  meridional  mass 
transport  after  60  years  integration  indicates  that  the  compatibility  consideration 
meridional  mass  transport  is  larger  than  that  of  incompatibility  consideration  in  all 
latitudes  of  the  northern  hemisphere  and  smaller  in  all  latitudes  of  southern  hemisphere. 
This  is  consistent  with  figure  4,  which  shows  an  extra  northwest  (southwest)  Ekman 
transport  in  northern  (southern)  hemisphere  when  the  compatible  forcing  is  used. 

Therefore,  we  see  the  difference  of  meridional  section  profile  of  v  velocity  in 
figure  24,  26,  and  28,  the  v  velocity,  which  is  due  to  compatible  forcing  correction,  is 
always  stronger  in  the  northern-hemisphere  than  that  due  to  incompatible  forcing 
influence.  In  the  southern-hemisphere  the  resulting  v  velocity  which  is  integrated  by 
compatible  forcing  is  found  to  be  weak.  The  errors  by  incompatible  forcing  is  most 
obvious  in  the  upper  sea  water.    The  errors  is  nearly  25  % . 

6.      Meridional  section  profile  of  w  velocity 

The  meridional  section  profile  of  w  velocity  at  160°E  after  60  years  is  shown 
in  figure  29,  and  the  difference  is  provided  in  figure  30.  The  meridional  section  profile 
of  w  velocity  at  160°W  after  60  years  is  shown  in  figure  31,  and  the  difference  is 
provided  in  figure  32.  The  meridional  section  profile  of  w  velocity  at  120°W  after  60 
years  is  shown  in  figure  33,  and  the  difference  is  provided  in  figure  34. 

From  the  meridional  section  profiles  of  w  velocity  in  figure  29,  31,  and  33, 
the  w  velocity  has  large  variation  in  longitude  (i.e.,  section  to  section)  indicating  the 
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effect  of  transients  (eddies).  The  w  velocity  profile  of  deep  water  is  more  complex  than 
that  of  upper  sea  water.  From  the  difference  of  meridional  section  profile  of  w  velocity 
in  figure  30,  32,  and  34,  the  w  velocity,  which  is  due  to  incompatible  forcing  influences, 
has  more  errors  near  the  equator  than  that  of  high  latitude.    The  errors  is  about  16%. 

7.      Difference  at  some  random  points 

a.  0  at  160°E,  0°N  on  level  2 

The  difference  in  6  for  depth  of  76.15  meters  (level  2)  at  160°E,  0°N 
is  shown  in  figure  35.  The  error  for  6  in  the  first  three  years  is  about  6%  due  to 
incompatible  forcing  influences.  This  error  increases  to  1%  after  about  10  years 
integration. 

An  extra  northwest  (southwest)  Ekman  transport  appears  in  the  northern 
(southern)  hemisphere  when  the  compatible  forcing  is  used.  Thus  an  additional 
horizontal  divergence  occurs  near  the  equator.  Causing  an  extra  amount  of  cold  deeper 
sea  water  moves  upward  toward  the  surface  making  the  temperature  lower  near  the 
equator  in  the  compatible  case.  Therefore,  at  depth  of  76.15  meters  (level  2)  at  160°E, 
0°N  the  d  by  the  compatible  forcing  in  the  whole  integration  time  of  60  years  is  always 
lower  than  that  due  to  incompatible  forcing  influences. 

b.  V  velocity  at  160°E,  5°S  on  level  2 

The  difference  in  v  velocity  for  depth  of  76. 15  meters  (level  2)  at  160°E, 
5°S  is  provided  in  figure  36.    The  error  in  v  velocity  in  the  first  three  years  is  about 
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26%  due  to  incompatible  forcing  influences.  This  error  increases  to  about  28%  after 
about  8  years  integration. 

From  the  meridional  section  profiles  of  v  velocity  in  figure  23,  25,  and 
27,  the  v  velocity  is  positive  from  the  equator  to  about  25  °N  where  the  velocity  changes 
direction  (opposite  sign)  toward  high  latitude.  The  v  velocity  in  the  southern-hemisphere 
is  compared  with  the  relation  in  the  northern-hemisphere. 

From  figure  14  and  16,  the  latitudinal  distribution  of  meridional  mass 
transport  after  60  years  integration  indicates  that  the  compatibility  consideration 
meridional  mass  transport  is  larger  than  that  of  incompatibility  consideration  in  all 
latitudes  of  the  northern  hemisphere  and  smaller  in  all  latitudes  of  southern  hemisphere. 
This  is  consistent  with  figure  4,  which  shows  an  extra  northwest  (southwest)  Ekman 
transport  in  northern  (southern)  hemisphere  when  the  compatible  forcing  is  used. 

Therefore,  at  depth  of  76.15  meters  (level  2)  at  160°E,  5°S  the  v 
velocity  by  the  compatible  forcing  in  the  whole  integration  time  of  60  years  is  always 
weaker  than  that  due  to  incompatible  forcing  influences. 

c.      V  velocity  at  160° E,  10° S  on  level  2 

The  difference  in  v  velocity  for  depth  of  76. 15  meters  (level  2  at  160°E, 
10  °S  can  be  seen  in  figure  37.  The  error  in  v  velocity  in  the  first  three  years  is  about 
23%  which  is  due  to  incompatible  forcing  influences.  This  increases  the  error  to  28% 
after  about  8  years  integration.  At  depth  of  76.15  meters  (level  2)  at  160°E,  10°S  the 
v  velocity  by  the  compatible  forcing  in  the  whole  integration  time  of  60  years  is  always 
weaker  than  that  due  to  incompatible  forcing  influences. 
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d.      W  velocity  at  160° E,  10° TV  on  level  5 

The  difference  in  w  velocity  for  depth  of  385.64  meters  (level  5)  at 
160°E,  10°N  is  illustrated  in  figure  38.  The  error  of  w  velocity  in  the  first  three  years 
is  about  16%  due  to  incompatible  forcing  influences.  This  error  decreases  to  10%  after 
about  8  years  integration.  At  depth  of  385.64  meters  (level  5)  at  160°E,  10°N  the  w 
velocity  by  the  compatible  forcing  is  always  stronger  than  that  due  to  incompatible 
forcing  influences  after  6  years  integration. 
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Figure  la.  The   surface   forcing   fields   which   were   used   by   Marotzke   and 

Willebrand  (1991). 
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Figure  lb. 


Zonal  mean  values  of  observed  sea  surface  temperature  (°C),  used  by 
Marotzke  and  Willebrand  (1991). 
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Figure  lc. 


Zonal  mean  values  of  observed  salinity  (psu),  used  by  Marotzke  and 
Willebrand  (1991). 
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Figure  2. 


Surface  zonal  wind  stresses  (dynes/cm2)  as  functions  of  latitude, 
(a)  without  thermally  driven  surface  winds  (dash  curve),  (b)  with 
thermally  driven  surface  winds  (solid  curve). 
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Figure  3. 


Surface  meridional  wind  stresses  (dynes/cm2)  as  functions  of 
latitude,  (a)  without  thermally  driven  surface  winds  (dashed 
curve),  (b)  with  thermally  driven  surface  winds  (dashed  curve). 
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Figure  4. 


The  Ekman  current  change  near  25  °N  and  25  °S  due  to 
compatible  forcing.  The  solid  arrows  show  the  change  in  the 
wind  stress  (Figs.  2  and  3)  and  the  double  arrows  show  the 
resulting  change  in  the  surface  Ekman  current. 
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Figure  5. 


Total  northward  transport  of  heat  after  3,   6,   9,    12  years 
integration,  (x  1015  Watts). 
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Figure  6. 


Difference  of  total  northward  transport  of  heat  after  3,  6,  9,  12 
years  integration,  (x  1015  Watts). 
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figure  7. 


Total  northward  transport  of  heat  after  51,  54,  57,  60  years 
integration,  (x  1015  Watts). 
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Figure  8. 


Difference  of  total  northward  transport  of  heat  after  51,  54,  57, 
60  years  integration  (x  1015  Watts). 
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figure  9. 


Total  northward  transport  of  salt  after  3,   6,   9,    12  years 
integration,  (x  1010  cm3/ sec). 
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Figure  10.  Difference  of  total  northward  transport  of  salt  after  3,  6,  9,  12  years 

integration,  (x  1010  cm3/ sec). 
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figure  11. 


Total   northward   transport   of  salt   after  51,    54,    57,    60 
integration,  (x  1010  cm3/sec). 
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Figure  12.  Difference  of  total  northward  transport  of  salt  after  51,  54,  57,  60 

years  integration,  (x  1010  crnVsec). 
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Figure  13.  Meridional  mass  transport  after  30  years  integration,    (in  Sverdrup). 
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Figure  14.  Difference  of  meridional  mass  transport  after  30  years  integration,  (in 

Sverdrup). 
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Figure  15. 


Meridional  mass  transport  after  60  years  integration,  (in  Sverdrup). 
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Figure  16.  Difference  of  meridional  mass  transport  after  60  years  integration,  (in 

Sverdrup). 
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Figure  17 


Distribution  of  temperature  (°C)  at  160°E  meridional  section  after  60 
years  integration. 


Figure  18.  Distribution  of  temperature  difference  (°C)  at   160°E  meridional 

section  after  60  years  integration  between  two  different  kinds  of 
surface  forcing. 
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Figure  19. 


Distribution  of  temperature  (°C)  at  160°W  meridional  section  after  60 
years  integration. 


Figure  20.  Distribution  of  temperature  difference  (°C)  at  160°W  meridional 

section  after  60  years  integration  between  two  different  kinds  of 
surface  forcing. 
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Figure  21.  Distribution  of  temperature  (°C)  at  120°W  meridional  section  after  60 

years  integration. 
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Figure  22.  Distribution  of  temperature  difference  (°C)  at  120°W  meridional 

section  after  60  years  integration  between  two  different  kinds  of 
surface  forcing. 
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Figure  23.  Distribution  of  v  (cm/sec)  at  160°E  meridional  section  after  60  years 

integration. 


Figure  24.  Distribution  of  v  difference  (cm/sec)  at  160°E  meridional  section  after 

60  years  integration  between  two  different  kinds  of  surface  forcing. 
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Figure  25.  Distribution  of  v  (cm/sec)  at  160°W  meridional  section  after  60  years 

integration. 
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Figure  26.  Distribution  of  v  difference  (cm/sec)  at  160°W  meridional  section 

after  60  years  integration  between  two  different  kinds  of  surface 
forcing. 
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Figure  27.  Distribution  of  v  (cm/sec)  at  120°W  meridional  section  after  60  years 

integration. 
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Figure  28.  Distribution  of  v  difference  (cm/sec)  at  120°W  meridional  section 

after  60  years  integration  between  two  different  kinds  of  surface 
forcing. 
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Figure  29. 


Distribution  of  w  ( x  10"3  cm/sec)  at  160°E  meridional  section  after  60 
years  integration. 


Figure  30.  Distribution  of  w  difference  (x   10"3  cm/sec)  at  160°E  meridional 

section  after  60  years  integration  between  two  different  kinds  of 
surface  forcing. 
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figure  31 


Distribution  of  w  (x  10~3  cm/ sec)  at  160°W  meridional  section  after 
60  years  integration. 


Figure  32 


Distribution  of  w  difference  (x  10"3  cm/sec)  at  160°W  meridiona 
section  after  60  years  integration  between  two  different  kinds  of 
surface  forcing. 
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Figure  33.  Distribution  of  w  (x  10"3  cm/sec)  at  120°W  meridional  section  after 

60  years  integration. 


Figure  34. 


Distribution  of  w  difference  (x   10"3  cm/sec)  at  120°W  meridiona 
section  after  60  years  integration  between  two  different  kinds  of 
surface  forcing. 
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Figure  35.  Difference  in  0  (°C)  for  depth  of  76.15  meters  (level  2)  at  160°E, 

0°N. 
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Figure  36.  Difference  in  v  (cm/sec)  for  depth  of  76. 15  meters  (level  2)  at  160°E, 

5°S. 
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Figure  37.  Difference  in  v  (cm/sec)  for  depth  of  76. 15  meters  (level  2)  at  160°E, 

10°S. 
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Figure  38.  Difference  in  w  (x  10"3  m/day)  for  depth  of  385.64  meters  (level  5) 

at  160°E,  10°N. 
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V.     CONCLUSIONS 

The  surface  wind  and  surface  thermal  conditions  (e.g. ,  SSAT,  net  surface  heat  flux) 
are  important  components  of  the  atmospheric  system  and  are  related  each  other. 
However,  the  atmospheric  parameters  representing  mechanical  and  buoyancy  forcing, 
such  as  surface  wind  and  SST  fields  (or  the  net  surface  heat  flux)  are  often  prescribed 
independently  in  current  ocean  models,  especially  idealized  or  "process"  models.  The 
practice  of  prescribing  wind  stress  and  buoyancy  forcing  independently  in  ocean 
numerical  models  is  an  unrealistic  assumption. 

A  large  difference  between  the  ocean  circulation  forced  by  compatible  and 
incompatible  surface  forcing  (i.e.,  surface  wind  stress  and  buoyancy  forcing)  has  been 
found  by  comparing  the  results  from  the  GFDL  MOM  model.  In  particular,  the 
compatible  forcing  derived  from  the  (incompatible)  forcings  field  of  Marotzke  and 
Willebrand  (1991),  produces  15-20%  stronger  meridional  circulation  cells  in  the  low 
latitudes  of  both  hemispheres  of  the  MOM  model.  This  difference  can  be  attributed  to 
the  incompatible  SSAT  (or  surface  heat  flux)  and  surface  wind  fields.  Since  most  ocean 
numerical  models  are  forced  models,  a  compatible  forcing  boundary  condition  is  an 
importance  consideration  before  running  any  numerical  ocean  model. 


59 


LIST  OF  REFERENCES 

1.  Chu,  P.C.,  1986:  An  Instability  Theory  of  Ice-air  Interaction  for  the  Migration  of 
the  Marginal  Ice  Zone.    Geophys.  J.  R.  Astr.  Soc,  86,  863-883. 

2.  Chu,  P.C.,  1989:  Relationship  between  Thermally  Forced  Surface  Wind  and  Sea 
Surface  Temperature  Gradient.    Pure  Appl.  Geophys.,  130,  31-45. 

3.  Chu,  P.C.,  1992:  Incompatible  Wind  and  Buoyancy  Forcing  in  Ocean  Numerical 
Models.   Adv.  Atmos.  Sci.,  (in  press). 

4.  Bryan,  K.,  and  M.  D.  Cox,  1967:  A  Numerical  Investigation  of  the  Oceanic 
General  Circulation.    Tellus,  19,  54-80. 

5.  Bryan,  K.,  1969:  A  Numerical  Method  for  the  Study  of  the  Circulation  of  the 
World  Ocean.   J.  Computat.  Phys.,  4,  347-376. 

6.  Cox,  M.D.,    1985:   An  Eddy  Resolving  Numerical  Model  of  the  Ventilated 
Thermocline.   J.  Phys.  Oceanogr.,  15,  1312-1324. 

7.  Cox,  M.D. ,  and  K.  Bryan,  1984:  A  Numerical  Model  for  Ventilated  Thermocline. 
J.  Phys.  Oceanogr.,  14,  674-687. 

8.  Davey,  M.K.,  1983:  A  Two-level  Model  of  a  Thermally  Forced  Ocean  Basin.  J. 
Phys.  Oceanogr.,  13,  169-190. 

9.  Haney,  R.L.,  1971:  Surface  Thermal  Boundary  Condition  for  Ocean  Circulation 
Models.   /.  Phys.  Oceanogr.,  1,  241-248. 

10.  Holton,  J.R.,  1972:  An  Introduction  to  Dynamic  Meteorology.    Academic  Press, 
New  York,  319p. 

1 1 .  Huang,  R.X. ,  1989:  Sensitivity  of  a  Multilayered  Ocean  General  Circulation  Model 
to  the  Surface  Thermal  Boundary  Condition.  J.  Geophys.  Res.,  94,  18011-18021. 

12.  Kuo,  H.L. ,  1973:  Planetary  Boundary  Layer  Flow  of  a  Stable  Atmosphere  over  the 
Globe.   J.  Atmos.  Sci.,  30,  53-65. 

13.  Marotzke,   J.,   and  J.   Willebrand,    1991:   Multiple  Equilibria  of  the  Global 
Thermohaline  Circulation.   J.  Phys.  Oceanogr.,  21,  1372-1385. 


60 


14.  Robinson,  A.R.,  D.E.  Harrison,  Y.  Mintz,  and  A.J.  Semtner,  Jr.,  1977:  Eddies 
and  the  General  Circulation  of  an  Idealized  Oceanic  gyre:  a  Wind  and  Thermally 
Driven  Primitive  Equation  Numerical  Experiment.  J.  Phys.  Oceanogr.,  7,  182- 
207. 

15.  Semtner,  A. J.,  Jr.,  and  Y.  Mintz,  1977:  Numerical  Simulation  of  the  Gulf  Stream 
and  Mid-ocean  Eddies.   J.  Phys.  Oceanogr. ,  7,  208-230. 

16.  Takano,  K.,  1975:  A  Numerical  Simulation  of  the  World  Ocean  Circulation: 
Preliminary  Results.  In:  Numerical  Models  of  Ocean  Circulation.  National 
Academy  of  Sciences,  Washington  D.C.,  121-129. 

17.  Willmont,  A.J.,  and  M.S.  Darby:  Nonlinear  Rossby  Waves  in  a  Thermodynamic 
Reduced  Gravity  Ocean.    Geophys.  Astrophys.  Fluid  Dyn.,  54,  189-227. 


61 


INITIAL  DISTRIBUTION  LIST 


1 .  Defense  Technical  Information  Center 
Cameron  Station 

Alexandria,  VA  22304-6145 

2.  Dudley  Knox  Library 
Code  0142 

Naval  Postgraduate  School 
Monterey,  CA  93943-5002 

3.  Chairman  Dr.  Robert  L.  Haney 
Code  MR/Hy 

Department  of  Meteorology 
Naval  Postgraduate  School 
Monterey,  CA  93943 

4.  Professor  Peter  C.  Chu 
Code  OC/Cu 

Department  of  Oceanography 
Naval  Postgraduate  School 
Monterey,  CA  93943 

5.  LT  Kuo,  Yu-Heng 

NO.  17-3,  Ping-Shan  Hsin  Tsun, 
Tsoying,    Kaohsiung, 
Taiwan,  R.  O.  C. 

6.  Naval  Academy  Library 
P.O.  Box  90175 
Tsoying,  Kaohsiung, 
Taiwan,  R.  O.  C. 


62 


DUDLEY  KNOX  LIBtv 

NAVAL  POSTGRADU-.        CHOOl 

MONTEREY  CA  9394^,01 


GAYLORD  S 


